First we load up the packages we’ll be working with
# remotes::install_github("mathesong/bloodstream")
# remotes::install_github("mathesong/kinfitr")
library(tidyverse)
── Attaching core tidyverse packages ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(bloodstream)
library(kinfitr)
library(job)
download.file("https://www.dropbox.com/scl/fi/ax1t3lfop3pp3h4cir0m7/ds004869.zip?rlkey=6z7e9wopud9ngu26ekbd3kfks&dl=1",
destfile = "ds004869.zip")
trying URL 'https://www.dropbox.com/scl/fi/ax1t3lfop3pp3h4cir0m7/ds004869.zip?rlkey=6z7e9wopud9ngu26ekbd3kfks&dl=1'
Content type 'application/binary' length 29402005 bytes (28.0 MB)
==================================================
downloaded 28.0 MB
unzip(zipfile = "ds004869.zip")
job({
bloodstream("ds004869/")
})
Alternatively, if you want to use a configuration file
bloodstream("ds004869/", configpath = "ds004869/code/config_tutorial.json")
We want the arterial input function data, which are called “input”.
bloodstream_data <- bloodstream_import_inputfunctions("ds004869/derivatives/bloodstreamtutorial/") %>%
select(-measurement)
head(bloodstream_data, n = 6)
bloodstream_data %>%
slice(1:6) %>%
pull(input) %>%
map(plot)
[[1]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[2]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[3]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[4]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[5]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
[[6]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_line()`).
petprep_data <- kinfitr::bids_parse_files("ds004869/derivatives/petprep_extract_tacs/") %>%
unnest(filedata) %>%
filter(str_detect(path_absolute, "gtmseg")) %>%
mutate(
ses = case_when(
ses == "01" ~ "baseline",
ses == "02" ~ "blocked"
)
)
tacs <- petprep_data %>%
filter(measurement=="tacs") %>%
filter(is.na(pvc)) %>%
mutate(tacdata = map(path_absolute, ~read_delim(.x, delim="\t",
show_col_types = FALSE)))
morphdata <- petprep_data %>%
filter(measurement=="morph") %>%
mutate(morphdata = map(path_absolute, ~read_delim(.x, delim="\t",
show_col_types = FALSE)))
Here, we would usually combine TACs according to which combined regions we were interested in, and use volume-weighted averaging. In this case, we’ll just skip that step and be a little bit naughty and focus on one region: the left putamen.
Now we want to combine regions by volume-weighted averages of the contitutent regions. We’ll create a frontal cortex region, a striatal region and a hippocampus-amygdala region.
First, let’s define the regions
frontal_regions <- selected_tacs$morphdata[[1]] %>%
filter(str_detect(name, "frontal")) %>%
pull(name)
frontal_regions
[1] "ctx-lh-caudalmiddlefrontal" "ctx-lh-lateralorbitofrontal" "ctx-lh-medialorbitofrontal" "ctx-lh-rostralmiddlefrontal"
[5] "ctx-lh-superiorfrontal" "ctx-lh-frontalpole" "ctx-rh-caudalmiddlefrontal" "ctx-rh-lateralorbitofrontal"
[9] "ctx-rh-medialorbitofrontal" "ctx-rh-rostralmiddlefrontal" "ctx-rh-superiorfrontal" "ctx-rh-frontalpole"
striatal_regions <- selected_tacs$morphdata[[1]] %>%
filter(str_detect(name, "Putamen|Accumbens|Caudate")) %>%
pull(name)
striatal_regions
[1] "Left-Caudate" "Left-Putamen" "Left-Accumbens-area" "Right-Caudate" "Right-Putamen" "Right-Accumbens-area"
hipamg_regions <- selected_tacs$morphdata[[1]] %>%
filter(str_detect(name, "Hippocampus|Amygdala")) %>%
pull(name)
hipamg_regions
[1] "Left-Hippocampus" "Left-Amygdala" "Right-Hippocampus" "Right-Amygdala"
Now we combine the TACs and the region sizes
selected_tacs <- select(tacs, c(ses:rec, tacdata)) %>%
inner_join(select(morphdata, c(ses:rec, morphdata))) %>%
group_by(sub, ses) %>%
mutate(tacdata = map(tacdata, ~pivot_longer(.x,
cols = `Left-Cerebral-White-Matter`:`ctx-rh-insula`,
names_to = "name", values_to = "TAC"))) %>%
mutate(tacdata = map2(tacdata, morphdata, ~inner_join(.x, .y, by="name")))
Joining with `by = join_by(ses, sub, task, trc, acq, run, rec)`
Then we perform volume-weighted averaging
out <- do.call(out, inner_join)
Error in do.call(out, inner_join) : second argument must be a list
modeldata <- selected_tacs %>%
select(ses:rec, tacdata = selected_tacdata) %>%
inner_join(bloodstream_data)
Joining with `by = join_by(ses, sub, task, trc, acq, run, rec)`
The TAC data are in Bq/mL, and the bloodstream data are in kBq/mL. So we correct this.
modeldata <- modeldata %>%
mutate(tacdata = map(tacdata, ~.x %>%
mutate(across(.cols = Frontal:HippAmg,
~unit_convert(.x, from_units = "Bq", to_units = "kBq")))))
modeldata <- modeldata %>%
mutate(tacdata = map(tacdata, ~.x %>%
mutate(frame_start = frame_start / 60,
frame_end = frame_end / 60) %>%
mutate(frame_dur = frame_end - frame_start,
frame_mid = frame_start + 0.5*frame_dur)))
modeldata <- modeldata %>%
mutate(tacdata = map(tacdata, ~.x %>%
mutate(meanTAC = rowMeans( .x %>% select(Frontal:HippAmg) )) %>%
mutate(weights = weights_create(t_start = frame_start,
t_end = frame_end,
tac = meanTAC))))
fit_2tc <- twotcm(
t_tac = modeldata$tacdata[[1]]$frame_mid,
tac = modeldata$tacdata[[1]]$Frontal,
input = modeldata$input[[1]],
weights = modeldata$tacdata[[1]]$weights,
vB = 0.05, multstart_iter = 5)
fit_2tc$par
plot(fit_2tc)
Here we’ll use a linearised model because they fit more quickly, in this case Logan. But most linearised models require a t* value to operate.
Let’s choose an appropriate t* value
modeldata %>%
filter(ses=="baseline") %>%
slice(1:5) %>%
mutate(tstarplot = map2(tacdata, input,
~Logan_tstar(
t_tac = .x$frame_mid,
lowroi = .x$HippAmg,
medroi = .x$Frontal,
highroi = .x$Striatum,
input = .y,
vB = 0.05)
)) %>%
pull(tstarplot)
Warning: There were 30 warnings in `mutate()`.
The first warning was:
ℹ In argument: `tstarplot = map2(...)`.
ℹ In group 1: `sub = "01"` and `ses = "baseline"`.
Caused by warning:
! Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 29 remaining warnings.
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
Ok, let’s use 13 frames.
Let’s focus on the frontal cortex
modeldata <- modeldata %>%
group_by(sub, ses) %>%
mutate(Logan_fit = map2(tacdata, input, ~Loganplot(t_tac = .x$frame_mid,
tac = .x$Frontal,
input= .y,
tstarIncludedFrames = 13,
weights = .x$weights,
vB = 0.05,
dur = .x$frame_dur)))
Let’s see a few fits
map(modeldata[1:6,]$Logan_fit, plot)
[[1]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
[[2]]
Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).
[[3]]
Warning: Removed 2 rows containing missing values or values outside the scale range (`geom_point()`).
[[4]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
[[5]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
[[6]]
Warning: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
Logan_outcomes <- modeldata %>%
select(sub, ses, Logan_fit) %>%
mutate(Vt = map_dbl(Logan_fit, c("par", "Vt"))) %>%
select(-Logan_fit)
ggplot(Logan_outcomes, aes(x=ses, y=Vt, colour=sub, group=sub)) +
geom_point(size=3, colour="black") +
geom_point(size=2.5) +
geom_line() +
scale_y_log10()